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by 

Juan G. Roederer 


SUMMARY 



A computer code is described which can be used for a two dimensional 
mapping of directional fluxes of energetic particles trapped in the outer mag- 
netosphere. G. Mead's model of the magnetospheric field is used. The main 
function of the code is to find, on a fixed reference meridian, the field line, 
equatorial pitch angle and other parameters of the shell generated by the par- 
ticles in their longitudinal drift, provided they are stably trapped. 


A series of curves is given, which can be used for a first rough, graphical 
attempt of flux mapping. 

The code is described in more detail in the Appendix. 
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A PROPOSAL FOR TRAPPED PARTICLE FLUX MAPPING 
IN THE OUTER MAGNETOSPHERE 


I. INTRODUCTION 

Several years ago, a considerable improvement in the study of the inner 
radiation belt was achieved by the introduction of a convenient, two dimensional 
mapping of trapped particles fluxes. 1 Two main features of the geomagnetic 
field in the inner magnetosphere made this possible: 

1. Time variations of the field within 2-3 earth radii are small enough, so 
that a world-wide adoption of a static field model is possible. 

2. The azimuthal asymmetry of the field is small enough, so that particle 
shells can be considered degenerate: 2 particles mirroring on a given 
field line populate identical shells, irrespective of their equatorial pitch 
angles, or their mirror points. Therefore, in a steady state, not only the 
directional flux j^, but also the omnidirectional flux J Q of trapped par- 
ticles is independent of longitude, being only a function of the mirror 
point field B m and the local L-value (and the energy): J Q = J Q (B m ,L,E) 

In the outer magnetosphere, none of these two characteristic features is 
true. The field configuration is strongly time dependent, driven by variations of 
the solar wind, by the ring current and by the diurnal wabbling of the earth's 
dipole. It further shows a strong day-night asymmetry, which causes consider- 
able shell splitting 3 . The fact, however, that some general features of the field 
configuration still seem to be preserved throughout all minor time variations 
may suffice to justify the use of a static model in a first attempt to organize 
particle flux data, measured at different positions in the outer magnetosphere 
during geomagnetic ally quiet epochs. In any case, such a procedure would be 
far better than using dipole quantities such as L-values, invariant latitudes, etc. 

Shell splitting does not allow a two dimensional description of omnidirectional 
particle fluxes. However, directional fluxes of trapped particles still can be 
mapped in a two dimensional way, independent of the local time at the point of 
measurement, if one makes use of a model for the magneto spheric field and 
computes theoretical shell configurations. Consider such a model, and take a 
directional flux measurement made at a point of geomagnetic coordinates R p , 

\ p , cp , of particles incident from a given direction (mean direction of viewing 
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of the detector). With the use of the model, one can determine the particles' 
pitch angle a p at the point of measurement, their mirror points and mirror field 
B , their I-value 

I =| VI -B(s)/E^ ds, 

* m 

equatorial pitch angle a e and radial distance to the equatorial crossing R e (in- 
tersection of the field line with the equatorial plane) . One now takes a fixed ref- 
erence meridian (for instance, the noon meridian), and computes with the aid of 
the model (and the conservation theorems of the adiabatic invariants) the equa- 
torial pitch angle a Q and the equatorial crossing R 0 of these particles, when 
they drift through this reference meridian, provided, of course, that they are 
stably trapped. The directional flux of the given class of particles at this ref- 
erence meridian will be the same as at the point of measurement, according to 
Liouville's theorem. This procedure leads to an analog to B-L plots of direc- 
tional fluxes, in which L is now replaced by the equatorial crossing R 0 in the 
refer'ence meridian. 

Magnetospheric models are not yet ripe for world-wide adoption; neverthe- 
less it seemed useful to propose here a computer code, which for a given model, 
allows two-dimensional mappings of directional fluxes of trapped particles, 
measured anywhere in the magnetosphere, and giving proper diagnostics if the 
particles were on open field lines, or if they were unable to reach the reference 
meridian in their longitudinal drift (pseudo-trapped particles *) . 

We want to emphasize again the limitations of such a procedure. However, 
we suggest that the adoption of such a "recipe" for flux mapping would be far 
more legitimate than the forced use of dipole relations in a region of the mag- 
netosphere which has little in common with a dipole field. 


H. THE COMPUTER CODE "ORGAN” 

This computer code, consists of two parts: the first one finds field line and 
shell parameters at the point of measurement for the particles in question; the 
second one finds the field line on which the particles will drift through the refer- 
ence meridian (provided they are stably trapped) . 

The first part performs field line tracings to the conjugate point (in order 
to compute I), to the equator (in order to find a e and R e ) and to the intersection 
with the earth's surface. The second part locates the point of prefixed I, B m 
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values on the reference meridian, using a fast iteration process (subroutine 
SEARCH), and then finds a 0 , R 0 and the intersection of the corresponding field 
line with the earth's surface. 

The magnetic field configuration adopted in the present code is that given 
by Mead 5 (subroutine MODMAG). This model has a severe limitation in that it 
takes the dipole axis perpendicular to the ecliptic. However, the problem of the 
interaction of a steady flow of plasma with a tilted dipole has not yet been 
solved. 

There are four parameters which have to be defined in this model: the 
minimum distance from the center of the earth to the magnetopause, in the 
noon meridian (called FRONT in the code) ; the distances from the center of the 
earth to the close and to the far limits of the neutral sheet in the midnight meri- 
dian (called SHEET1 and SHEET2, respectively), and the field intensity near 
the neutral sheet (called FSHEET). Suggested values 6 for these parameters 
are FRONT = 10 earth radii, SHEET 1 = 8 earth radii,* SHEET2 = 200 earth 
radii and FSHEET =15 gammas. Time variations can be simulated by adequate 
changes of these parameters. 3 The general field configuration is mainly de- 
termined by the value of the parameter FRONT; the trapping boundary and the 
limits of the pseudo-trapping regions, 3 however, are strongly influenced by the 
tail parameters. 

Electric fields are not taken into account in this paper. This imposes the 
restriction of what follows, to energetic particles only. 

A description of the complete code is given in the Appendix . Input to the 
program are the geomagnetic coordinates of the point of measurement, and the 
local pitch angle (the code can easily be modified to read in a direction in space, 
computing the local pitch angle for the model in question) . Output are the geo- 
magnetic coordinates of the mirror point, equatorial crossing and intersection 
with the earth's surface of the field line going through the point of observation, 
and the same quantities at the reference meridian, as well as the mirror point 
field intensity (in gauss) and I (in earth radii). Diagnostics are given if particles 
are initially on an open field line, if they are trapped in secondary, off-equator 
field minima (called "minimum-B pockets" in the code), if they precipitate into 
the ionosphere, or if they are pseudo-trapped, i.e., unable to complete a 180° 
drift around the earth. Computer time is in general of the order of a few sec- 
onds per case, on an IBM 7094. However, a case of pseudo-trapping may take 
as long as 15-20 seconds (subroutine SEARCH must make sure that there exists 


*This does not necessarily mean that the neutral sheet actually starts at 8 earth radii; field lines 
are still closed at 15—18 R e in the equatorial plane in the tail, in this case. 
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no point on the reference meridian with the wanted pair of I, B m values). An 
option is given, which repeats all computations, now taking the midnight meridian 
as reference (this is necessary if one wants to know whether the particles are 
stably trapped) . 


III. GRAPHICAL RESULTS 

Although the computer code presented here is designed to be applied in- 
dividually for each point of measurement and each direction in space (i.e., 
pitch angle) , we now present a series of graphs which may be used for a first, 
very rough attempt to organize directional flux data in the outer magnetosphere. 

All graphs shown in Figure 1 display equatorial pitch angle a 0 vs. equatorial 
crossing R Q in the reference meridian (noon); each graph corresponds to a fixed 
local time at the point of measurement. There are three groups of curves in 
each graph, corresponding to three different values of initial pitch angles at the 
point of measurement (full lines: a p = 90°; broken lines: a p = 60°; dotted lines: 
a = 30°). The curves shown are of constant latitude and of constant radial dis- 

P 

tance of the point of measurement. Points of intersection of equal latitude and 
radial distance in each group are joined by thin lines; these thin lines therefore 
represent the relation of a Q and R Q for particles passing through the same 
point, situated at a given local time, with different pitch angles. There are 
regions in these graphs which are "inaccessible," i.e., for which no stably 
trapped particles exist. 

An example is now given on how to use these graphs. Suppose that a 
directional particle flux measurement was made at a point at 7 earth radii and 
20° geomagnetic latitude, at 0400 hours local time. The local average pitch 
angle was 60° . In order to find the equatorial pitch angle a Q of these particles 
at the reference meridian, and their equatorial crossing R q , one goes to the 
0400 LT graph, sorts out the 60° pitch angle set (broken lines), and finds the 
intersection of the a = 20°, R = 7 earth radii curves. For intermediate 

P P 

values, an interpolation scheme can be used. 

In order to facilitate conversion of these a 0 - R 0 plots into B m - R 0 plots 
(equivalent to the B-L coordinates), we present in Figure 2 the relation of mir- 
ror point field vs. equatorial pitch angle a Q , for different equatorial crossing 
R 0 in the noon meridian. 

We leave it up to the reader to find out the main features revealed by these 
graphs, in particular, to inspect how and where the azimuthal asymmetry causes 
the greatest deformation of particle shells. For a more physical discussion of 
all these effects, see Ref. 3. 
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APPENDIX 


Main Program ORGAN 

The main features of this program were explained in the text (Section II) . 
Figure 3 shows a simplified flow chart. More information can be obtained from 
the comment cards in the listing given at the end of the Appendix, in particular 
regarding input. 

Diagnostics are, in general, self-explanatory. "Minimum-B pocket" is the 
name of a region off the equatorial plane, in which the absolute field intensity 
passes through a secondary minimum along a field line. Two of such regions at 
positive and negative latitudes in the front side, near the boundary, are always 
present in a magnetospheric model; they also can show up in the tail, depending 
on the particular value given to the tail parameters. Whether they really occur 
in the magnetosphere is still an experimentally unsolved question. 


Subroutine SEARCH 


This subroutine consists of an iteration process which finds a point P with 
prefixed B and I values at a specified longitude. In addition, this subroutine 
gives detailed information about the field line which passes through the point P. 
Input for the subroutine are the geomagnetic coordinates of a starting point for 
the iteration process, the specified B and I values, and the prefixed tolerances 
for each of B and I. Output are the coordinates of the desired B-I point, and the 
coordinates of a series of points along the field line between the point P and its 
conjugate. 

There are two basic processes in this subroutine: one varies the radial 
distance searching along a constant latitude and longitude line for a point of 
given B; the other searches along a nearly constant- B line in a given meridian 
plane, for a point of prefixed I. The first of these processes is very fast; it only 
involves calling the subroutine MODMAG which computes the magnetic field at 
each point. The second process is much more complicated and involves the use 
of subroutine INVAR, which traces a field line through each point at which it is 
called and evaluates the integral I. The program was constructed in such a way 
as to minimize the number of times INVAR is called. 

Figure 4 shows a simplified flow chart for subroutine SEARCH. Geo- 
metrically, the procedure is as follows: starting at the input point, the program 
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iterates the radial distance to a point Q where the prefixed B-value (called SB) 
is reached within a prefixed, gross tolerance designated by ERRB. The I value is 
then determined for this point Q. If it is not equal to the prefixed I-value (called 
SI) within the tolerance designated by ERRI, the latitude is changed by a small, 
prefixed amount DV in the correct direction. At this new latitude, a point R is 
found with the prefixed B-value; I is then calculated at this point. Points Q and 
R lie., approximately, on a constant B line in the meridian plane. A linear ex- 
trapolation is then used to find the point on this line which would have the speci- 
fied value SI. This extrapolation is made as though the I-dependence were linear 
along the constant B line. This extrapolated point is now taken as a second ap- 
proximation and is used as the starting point of the next iteration cycle. If the 
new starting point is close enough to the points Q and R, the step size in latitude 
DV is reduced by one-tenth. The iteration is continued until a point is found 
whose B and I values lie in the interval SB ± ERRB and SI ± ERRI. When this is 
accomplished, all tolerances are reduced two orders of magnitude and the entire 
procedure is repeated once more until a point is found whose B and I values 
agree with the prefixed ones within the new, reduced tolerance intervals. 

If iterations exceed a given number of cycles (set = 20 in this code) , a cut- 
off is introduced and a diagnostic is printed, telling that the code is unable to 
find a point with the specified I,B values. It should be noted that under normal 
conditions, it never takes more than 2-3 cycles in each iteration. 

Notice that subroutine SEARCH can be used to determine B-I rings in the 
outer magnetosphere. 


Subroutines INVAR, START, LINES, INTEG and MODMAG 

The first four subroutines are slightly modified and implemented versions 
of Mcllwain's program for L-computations. 7 Their purpose is to trace a field 
line (START and LINES with MODMAG) and to determine the second invariant I 
(INTEG). Input for INVAR is a point P; output are the B and I values at P, and 
the coordinates of a sequence of points along the field line beginning at P and 
progressing to its conjugate. 

INVAR controls the rest of the above listed programs. When called from 
INVAR, subroutine START picks the first three points of a sequence along the 
field line, which are such that the second point is P and the B-values of the 
three are in decreasing order. LINES is a tracing routine which continues the 
sequence of points along the field line in the direction defined by START. This 
sequence progresses to the first point for which the B-value is again equal or 
exceeds that of P. The arc length between consecutive points in the sequence is, 
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Figure 


3_Flow chart for program ORGAN. 


9 













4 



Figure 3(continued)-Flow chart for program ORGAN. 
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Figure 4-Flow chart for subroutine SEARCH 
















in general, a linear function of the geocentric distance of the points, initially 
being directly controlled by the prefixed value of the parameter called ERR. 
During the tracing, the arc length keeps doubling until an error control (dependent 
on ERR) levels it off at a value which is usually eight or ten times greater than 
the initial size. The result is speed and accuracy at the expense of a nonconstant 
and sometimes unpredictable cell size. The main merit of Mcllwain's tracing 
routine clearly shows up when it is used for a magnetospheric model, where 
there are sudden, sharp bends of field lines near the boundary, and across the 
neutral sheet: in these regions, the automatic error control in LINES drastically 
reduces the cell size to ensure proper tracing. 

Subroutine INTEG was left intact from Mcllwain's code. Subroute MODMAG 
was written by G. Mead 8 and corresponds to his magnetospheric model . 5 This 
model is symmetric with respect to the equatorial plane, and field line tracings 
to determine I values could be performed from the point in question to the equator 
only. This time saving modification, however, was not introduced in the present 
code, in order to leave it flexible enough to be used with unsymmetric models. 


Subroutine EQUAT 

This subroutine commands field line tracing from a given point to the geo- 
magnetic equator, i.e., to the point with minimum B value. Usually, the starting 
point is one defined in a previous INVAR call and is already close to the equator. 
EQUAT is always used for a more accurate determination of the equatorial 
point and is called with a much smaller value of ERR. When called from EQUAT, 
the cell-doubling mechanism in LINES is by-passed to ensure highest precision. 


Subroutine BESECT 

This subroutine operates similarly to EQUAT, but traces a field line down- 
wards, i.e., in the direction opposite to the equatorial point, until a prefixed 
B-value is reached. 


Subroutine INSECT 

This subroutine traces a field line from a given point downwards to the 
intersection with a fixed altitude (the earth's surface). 

All operations mentioned in connection with subroutines EQUAT, BESECT 
and INSECT are actually executed in START, LINES and MODMAG. Information 
about the particular type of tracing is transferred by a parameter called MMM. 
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PROGRAM ORGAN 

COMMON FR0NT,SHEfcTl.SHbET2.FSHEET 

COMMON a<2rjn).VNl<2O0).VN2«200>»VN3(200 >*ARC(200 )»VNEAR(3)*VBO<3)# 
1 VSaVE<3)#80#BNEaW*JUP#MMM 
COMMON VC0NJ(3) , alto 

TOR ENQUIRIES about this CODE, write to wUAM g. roederer 
FACUHAD DE CIEN01AS EXACTAS, PERU 272 , 0UENQS AIRES ARGENTINA 

DIMENSION V ( 3 ) 

I F LESS ACCURACY BUT HIGHER SPEED IS WANTED, MULTIPLY ALL ERRORS 
BY A COMMON FACTOR (NOT EXCEEDING 40-50) 

ENR=o . 0001 
E R R B = 0 .004 
ERRI=t),006 
ERR1=£RR*G . 01 
ERR2 = bRR*40 . 

ERR3 = ERR*0 . 05 
ZERO = 0 . 

R A D = 5 7 ,2957795 
RADIN=1./RAD 
AL TO = 1 . 

FRONT IS DISTANCE IN EARTH RADII TO M AGNE TOSPHER I C BOUNDARY IN THE 
EARTH SUN DIRECTION (SUGGESTED VALUE ? 10 ) 

SHEET1 IS DISTANCE TO INNER EDGE OF NEUTRAL SHEET IN midnight 
MERIDIAN (SUGGESTED VALUE 3 8-10) 

SHEET 2 IS DISTANCE TO OUTER EDGE OF NEUTRAL SHEET IN MIDNIGHT 
MERIDIAN (SUGGESTED VALUE 3 200) 

FSHEET IS MAGNETIC FIELD INTENSITY IN GAMMA NEAR NEUTRAL SHEET 
(SUGGESTED VALUE 3 15 GAMMA) 

READ (5.17) FRONT, SHEEU, SHEET*, FSHEET 
WRITE (6,17) FRONT, SHEETl,SHEET2, FSHEET 

IF BOTH NOUN AND MIDNIGHT MERIDIANS ARE WANTED AS REFERENCE PLANES 
SET NREF = 1 , t F ONLY NOON MERIDIAN IS WANTED AS REFERENCE, SET 

nrff = o 

R t AD (5,18) NREF 

ALT = GEOCENTRIC DISTANCE IN EARTH RADII TO POINT OF OBSERVATION 
FLAT = GEOMAGNETIC LATITUDE O r SAME POINT 

FLONG = LOCAL TIME OF SAME POINT, IN DEGREES EAST OF MIDNIGHT 
(E.G. DAWN MERIDIAN HAS *90 DEG. , DUSK MERIDIAN -90 OR *27Q ) 

alpha = local pitch angle in degrees, code may be easily implem* 

ENTEO SUCH AS TO READ A DIRECTION IN A GEOMAGNETIC SYSTEM, 
COMPUTING THE THEORETICAL PITCH ANGLE ASSOCIATED TO THIS DIRECTION 
in place indicated in comment cards below 

read (5,17) ALT, flat, flong, alpha 
IF (ALT.LT. 0.5) GO TO 16 
WRITE (6,23> ALT, f- LAT, FLONG, Alpha 
REFLON = 180 . 

IP = 0 

I F (aBS(FlAT) .LT. l. ) F L A T = 0 .4 

V(1>=ALT 


A 1 
A 2 
A 3 
A 4 
A 5 
A 6 
A 7 
A 8 
A 9 
A 10 
A 11 
A 12 
A 13 
A 14 
A 15 
A 16 
A 17 
A 18 
A 19 
A 20 
A 21 
A 22 
A 23 
A 24 
A 25 
A 26 
A 27 
A 28 
A 29 
A 30 
A 31 
A 32 
A 33 
A 34 
A 35 
A 36 
A 37 
A 38 
A 39 
A 40 
A 41 
A 42 
A 43 
A 44 
A 45 
A 46 
A 47 
A 48 
A 49 
A 50 
A 51 
A 52 
A 53 
A 54 
A 55 
A 56 
A 57 
A 58 
A 59 
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V(2)*<90.-FLAT)*RADIN A 60 

V(3)=FL0NG*RADIN A 61 

' SIT*aBS(SIN(V(2>>) ' ' * 62 

A 63 

THE FOLLOWING PART FINDS THE MAGNETOSPHER IC GEOMETRIC PARAMETERS A 64 
FOR THE FIELD LINE THROUGH THE POINT OF OBSERVATION A 65 

A 66 

CALL MODMAG ( ALT # S I T , V ( 3 ) . BR , BT , BP, 8B » V < 2 ) > A 67 

IF (aLPHA.LT. 0.5) ALPHA*0.5 A 68 

SINAsSINt ALPHA*RADIN) A 69 

SSQA=SINA*SINA A 70 

BMIR=8B/SSQA A 71 

BOsBMIR A 72 

A 73 

BESECT FINDS "MIRROR POINT ON INITIAL FIELD LINE ' '*'"74 

A 75 

CALL BESECT (V<1>.V(2),V(3), BALT. BLAT.8L0NG. ERR) A 76 

SAlTsBALT A 77 

SLaTsBLAT A 78 

SLONG = BLONG _ * 79 

- IT' (SAL'T . GT'. ALD' GO TO 6 """ " ' * 81) 

A 81 

INVAR DETERMINES I VALUE A 82 

A 83 

CALL INVAR ( VBOd ) » VBO< 2 ) » VB0(3) , ERR» BMIR.F I ) A 84 

IF PRINTOUT OF ALL FIELD LINE POINTS BETWEEN MIRROR POINTS IS A 85 

“ WANTED. INSERT FOLLOWING CARDS “ ~ " * 86 " 

DO 500 J*2.JUP A 87 

OLAT=90.-VN2( J)*RAD A 88 

0L0NG=VN3( J)*RAD A 89 

IF (OLONG.GT.180. ) OLONGeOLONG-360 . A 90 

IF (OLONG.LT, -180. ) OLONG=OLONG*360. A 91 

50 0 WRITE ( 6.501 ) J.VNK J),OlAT,OlONG ' 92" 

501 FORMAT ( I5.3F15.3.E20.4) A 93 

IF ( VBO( 1 ) . GT . 1 , ) GO TO 2 A 94 

WRITE (6.20) A 95 

IPsl A 96 

VBO ( 1 ) = V ( 1 ) A 97 

VB0(2) = V(2) ~ " A 9G — 

VBO ( 3 ) =V ( 3 ) A 99 

IF (FI .GT.40. > GO TO 7 A 100 

DO 3 J=1,JUP A 101 

IF ( VNl ( J) . GT . 20 • ) GO TO 7 A 102 

_ CONTINUE _ _ A 103 

EOUAT FINDS EQUATORIAL POINT ON INITIAL FIELD LINE A 105 

A 106 

CALL fcQUAT (VNEaR(1>.VnEaR(2).VNEAR(3).E0B.EALT,ELAT.EL0NG,ERR1) A 107 

IF ( ABS(ELAT) .GT.5. ) GO TO 4 A 108 

RATIO=EQB/BMIR A 109 

- TGT = SORT( RATIO /(I .-RATIO)) “ A 118 

PITCHsATAN(TGT)*RAD A 111 

WRITE (6,24) EALT.ELONG.EQB, PITCH A 112 

GO TO 5 A 113 

4 WRITE (6.19) cALT.ELAT.ELONG A 114 

5 CALL INSECT (VBO(1).VBO(2).VBO<3),CALT,CLAT.CLONG,ERR3) A 115 

IF (IP.EQ.O) WRITE (6,25) S ALT , SLAT , SLONG , BMIR ,F r " ‘ " A 116 

WRITE (6.26) CLAT.CLONG A 117 

IF (IP.EQ.1) GO TO 1 A 118 

KOUNT=l A 119 
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GO TO 8 

WRITE (6.31) SALT. SLAT, SLONG.0MIR 

CALL INSECT (VB0<1).VB0<2).V80<3),CALT,CLAT.CL0NG.ERR3) 

WRITE (6.26) CLAT.CLONG 
GO TO 1 

call INSECT (VBO<1).V0O(2>.VBO(3).CALT.CLAT.CLONG.ERR3> 

WRITE (6,22) 

IF (IP.EQ.O) WRITE (6.25) SALT.SLAT.SLONG.BMIR.fi 
write (6.26) CLAT.CLONG 
GO TO 1 

THE FOLLOWING PART FINDS THE MAGNETIC PARAMETERS ON THE NOON 

and/or the midnight meridian 

COORi*BALT 

C00R2*BLAT 

C00R3=REFL0N 

search looks for a PREFIXED B-L POINT 

'CAU'SEARCH (COOftl.COOR2,COOR3.BMIR,FI.BBB.FFI,ERR2.ERRB,ERRI.CHEC 
IK) 

IF <<CHECK.GT.0.5).OR. ((KOUNT.EQ.l). AND. <VNEAR(1 >.GT. 10.5))) GO TO 
1 12 

IF PRINTOUT or ALL FIELD line points between MIRROR POINTS IS 

wanted, INSERT FOLLOWING CARDS _ _ 

DO 502 J«2. JUP' 

OLAT=90.-VN2< J)*RAD 


A 120 
A 121 
A 122 
A 123 
A 124 
A 125 
A 126 
A 127 
A 128 
A 129 
A 130 
A 131 
A 132 
A 133 
V 134 
A 135 
A 136 
A 137 
A 138 
_A _139 
A 140 
A 141 
A 142 
A 143 
A 144 
A 145 
A 146 
A 147 


0L0NG»VN3( J)*RAD * 148 

IF (OLONG . GT , 180 • ) OLQNGsOLONG-360 . A 149 

IF (OLONG. LT. -180. ) OLONG*OLONG*360. A 150 

502 WR I TE ( 6. 501 ) J . VN1 ( J >. OL AT . OLONG A 151 

CALL EOUAT (VNEAR(l).VNEAR(2),VNEAR<3).RE0B.REALT,RELAT7R'EL0NG.ERfr A 152 

11) A 153 

CALL INSECT ( VSA VE ( 1 ) , VSA VE ( 2 ) . VSA VE < 3 ) , C ALT , CL AT . CLONG, ERR3 > A 154 

IF (aBS(RELAT) .GT.5. ) GO TO 10 A 155 

RAT I O S REOB/0M I R A 156 

TGT*SORT(RATIO/(l. -RATIO)) A 1 57 

RALPHA»ATAN(TGf)*RAD " ' A 138 

IF (K0UNT.EQ.2) GO TO 9 A 159 

WRITE (6.27) REALT,REOB,RALPHA.COOR1.COOR2.CLAT A 160 

GO TO 14 A 161 

9 WRITE (6.28) REALT.REQB.RALPHA.COORl.COOR2.CLAT A 162 

GO TO 14 A 163 

10 “WRITE (6,2l) REaLT'.RELAT, RElONG * 1*64 

IF (K0UNT.EQ.2) GO TO 11 A 165 

WRITE (6.27) ZERO. ZERO. ZERO, COORl.COOR2.CLAT A 166 

GO TO 14 A 167 

11 WRITE (6.28) ZERO. ZERO. ZERO, COORl.COOR2.CLAT A 168 

GO TO 14 A 169 

12" If (KOUNT.EO.T) GO TO 13 ~ "A" 170 

WRITE (6.29) A 171 

GO TO l A 172 

13 WRITE (6,30) A 173 

GO TO 15 A 174 

i! IF_(K0UNT.EQ 1 2)_ GO TO 1 A 175 

15 IF (nReF.EQ.O)' GO TO 1 " A 176 

K0UNT*2 A 177 


REFLONa 0 , A 178 

GO TO 8 A 179 
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16 CONTINUE A 130 

RETURN A 181 

C ' “ A 182 

17 FORMAT (4F10.3) A 183 

18 FORMAT (15) A 184 

19 FORMAT QX/37H MIN1MUM-B POCKET IN INITIAL LINE AT.5X.4HRE *,F7.3 A 185 

1,5X.4HLATs,F6.2.5X,5HL0NGs.F7,2> A 186 

20 format (ix/39h particles precipitate into atmosphere) a i87 

2 1 ~ format (1X/35H MINIHUM-B POCKET JN FINAL LINE AT , 7X , 4HRF a , F71 3, 5 A 188 

1X#4HLATs,F6.2»5X*5HL0NGs*F7,2) A 189 

22 FORMAT (1X/30H PARTICLES ON OPEN FIELD LIME/) A 190 

23 FORMAT (1X////10H NEW CASE,9X,4HRE a , F7 . 3 , 5X , 4HL AT= , F6 . 2 , 5X , 5HL0N A 191 

1G=,F7.2,5X,12HPITCH ANGLE®, F5.2///1X) A 192 

24 FORMAT (9H EQUATOR, 10X.4HRE * , F 7 . 3 , 5X , 5HL0NG* , F7 . 2 , 5X , 3H Ba,F7.6, A 193 

15XV6hPITCH*,F6.2) A 194' 

25 FORMAT (14H MIRROR PO I NT , 5X , 4HRE * . F 7 . 3 » 5X , 4HL AT* . F 6 • 2 . 5X , 5HL0NG* a 195 

1,F7.2.5X.3HBM*.F7.6,5X2HI»»F6.2) A 196 

26 FORMAT (27H INTERSECTION WITH SURF ACE, 8X , 4HLAT s , F6 . 2, 5X , SHLONGs , F A 197 

17.2) A 198 

27 FORMAT (1X//18H ON NOON MER I D I AN//8H EQUATOR, 12X. 2HR* , F7 . 3, 5X . 2HB A 199 

1=,F7.6,7X,6HPITcHs,F6.2/i 3H MIRROR POINT, 5X.4HRE a , F7 .-3", 5X , 4HL AT* , A 200~ 

2F6.2/13H INTERSECTION, 5X,4HLAT»,F6. 2) A 201 

28 FORMAT (1X/22H ON MIDNIGHT MERIDIAN//8H EQUATOR, 12X , 2HRs , F7 . 3, 5X , A 202 

12HB*,F7.6,7X,6HPITcHs,F 6.2/13H MIRROR POl NT , 5X , 4HRE a , F7 . 3, 5X , 4HL A A 203 

2T*,F6.2/13H INTERSECTION, 5X, 4HL ATs , F6 . 2) A 204 

29 FORMAT (1X/37H PSEUDO-TRAPPED# LEAVES THROUGH TAIL/) A 205 

30" FORMAT (1X/45H PSEUDO-TRAPPED* LEAVES' THROUGH MAGNETOPAUSE/ ) " ~ ~A 2C6 

31 FORMAT (1X/29H TRAPPED IN MINlMUM-B POCKET /14H MIRROR POINT, 5X, 4 A 207 

1HRE =,F7.3,5X,4HLATa,F6.2.5X,5HL0NGs.F7.2.5X,3HBM*,F7.6) A 208 

END A 209- 
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SUBROUTINE SEARCH ( alt, Fl AT,FLONG,SB f SI »BB»F I . ERR,ERRB,ERRI , RLOST ) 
COMMON FRONT, SHEETi.SHEET2#rSHEET 

COMMON B(2n0).VNl<200).VN2(200),VN3(200)»ARC(200),VNEAR(3),VBO(3). 
1 VSaVE(3),80.BNEaR. JUP.MMM 
COMMON VCON J ( 3 ) » A|_ TO 

subroutine defines field LINE going THROUGH POINT OF prefixed 
B AND SECOND INVARIANT I , AT A GIVEN LONGITUDE 

DIMENSION V( 3). V 1(3), V2<3> 

DV=0 . 02 
RLOST = 0 . 

MCHECKsO 

I CHECK® 0 

SERR=ERR 
SERRB=ERRB 
SERR I *ERR l 
V(1)=ALT 

IF <V<1).LT.0.5) GO TO 11 _____ 

V< 2>=< 90. -FLAT J/57. 2957795 
V(3)=FL0N3/57. 2957795 
UCLT=1.570« 

I CON= 1 

I L I T = 1 

DELV2=DV 

iCHtcKs ICHECK+1 

IF IICHECK.GT.20) GO TO 12 

SIT*aBS(SIN(V<2) ) ) 

CALL MODMAG ( V ( 1 ) , S I T, V < 3 ) $ 8R , BT , BP, BB, V ( 2 ) > 

FAC*i.-(SB-BB)/(3.*S8) 

IF (FAC.GT.1.5) F AC S 1 . 5 
IF (FAc.Lt . 0 .666) F AC 3 0 . 666 
V(1)=V(1)*FAC 

IF {( V(l) .GT.100 • > -OR* (V(l) .LT.O .5) ) GO TO 11 
VKD'VIl) 

V1(2)=V(2) 

MCHECK=MCHECK*1 

IF (MCHECK.GT.15) GO TO 13 

CALL MODMAG ( V ( 1 ) , S I T , V ( 3 ) . BR , BT , BP , BB , V ( 2 > ) 

IF (aBS((8B-SB)/SB).GT.ERRB) GO TO 3 
MCHECK=0 

IF (ILIT.NE.1) GO TO 7 
ILIT=2 

CALL INVAR ( V ( 1 ) . V ( 2 ) . V (3 ) . ERR * BB , F 1) 

IF (JUP.LT.O) GO 10 14 
V2(1)»V(1) 

V2(2)*V<2) 

B2 = BB 
Fl2sFl 

IF (ABS( (FI-SI ) /S I ) .LE.ERRI ) GO TO 8 * 

IF (AbS(V(2)-DCLT).LT,0.1) GO To 6 
SGNsS I GN ( 1 . , (FI-SI ) ) 

IF ( V < 2 ) .LT.DCLT) GO TO 4 

DElV2 3 -SGN*DELV2 

GO TO 5 

DELV2=SGN*DELV2 

V(2)=V(2)+DELV2 
GO TO 2 

V<2) = V(2)* DEL V2 


B 1 
B 2 
B 3 
B 4 
B 5 
B 6 
B 7 
B 8 
B 9 
B 10 
B 11 
B 12 
B 13 
B" 14 
B 15 
B 16 
B 17 
B 18 
8 19 

"B 20 
B 21 
B 22 
B 23 
B 24 
B 25 
B 26 
8 27 

B 28 
8 29 

B 30 
B 31 
B 32 
B 33 
B 34 
B 35 
B 36 
B 37 
B 38 
B 39 
B 40 
B 41 
B 42 
B 43 
8 44 

B 45 
B 46 
B 47 
B 48 
a 49 
B 50 
B 51 
B 52 
B 53 
B 54 
B 55 
B 56 
B 57 
B 58 
B 59 
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10 

11 

12 

13 

14 

15 

16 


C 

C 

C 

17 

18 

19 

20 


CALL INVAR (V(l>'V<2)#V(3)*ERR#8B,n> 

IF (JUP.LT.O) GO TO 14 

IF <<H-sn*<FI-Fl2>.LE.0.> GO T 0 2 

V(2)=V(2)-2.*DELV2 

GO TO 2 

B 1 - 0 B 

CALL INVAR (V(1)#V(2),V(3)#ERR#BB»FI ) 

IF (JUP.LT.O) GO TO 14 
IF UBS((FJ-SI>/S1).LE.ERRI> GO TO 8 
FACT=(SI-FI ) / ( F 1 2-F 1 ) 

if ( abs(fact) , gt • 3 . ) fact =3. *sign<i. # fact > 

VC DsVKDM V2(1>-V1(1))#FACT 
V(2)=V1(2>^(V2(2)-V1(2))*FACT 
ys AH I NIC ABS( V(2)-V1(2) ), A0S ( V ( 2 ) - V2 < 2 ) ) ) 

IF CY.'Gt.ABSCVl(2>-V2<2) >) GO TO 1 
DV = Y 
GO TO 1 

continue 

IF ( IC0N.EQ.2) GO TO 9 

IC0Ns2 

DV=DV*0 . 1 

ERR=ERR*0 . 025 

ERRB=ERR8*0 .05 

ERR I =ERR I *0 .04 

GO TO 1 

ALT*V ( 1 ) 

FLaT=90.-V(2)*57.2957795 
FL0NG=V<3)*57. 2957795 
IF (FLONG.GT.180. ) FLONG=FLONG- 360 . 

IF (FLONG.LT. (-180 .) ) FLONGaFLONG+360* 

DO 10 I »1 , 3 
VSaVEC I )*VC I ) 

GO TO 16 

WRITE (6*17) VCl>,FLAT,FLONG 
GO TO 15 

WRITE (6,18) ICH&CK 
GO TO 15 

WRI TE ( 6i 19 ) ICHECIOMCHECK 
GO TO 15 

WRITE (6,20) ICHECK,MCHECK#ALT#FLAT,FLONG 
JUP = 1 
RLOST = 1 • 

ERRsSERR 

ERRBsSERRB 

ERR I rSERRI — 

RETURN 


FORMAT ( 75X # 19H ALT OUT OF L I M I TS/75X , 3F 10 . 3 ) 

FORMAT (69X,51H SORRY, 0UT I CANNOT FlNO THAT DAMN “POINT “IN TCKECK 
1/69X , 2 1 10 ) 

FORMAT (69X,51H SORRY, BUT I CANNOT FIND THAT DAMN POINT IN MCHECK 
1/69X ,2110) 

FORMAT ( 75X , 40H SORRY, BUT POINT IS IN THAT DAMN P0CKET/55X,2IlQ>3 

1E15._5J _ 

END 


B 

B 

B 

B 

B 

B 

B 

B 

B 

B 

B 

B 

B 

B 

"B 

B 

B 

8 

B 

B 

“B“ 

B 

B 

B 

B 

B 


60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

Dtr 

81 

82 

83 

84 

85 


B — 66“ 
B 87 


B 

B 

B 

B 

B 

B 

B 

B 

B 

B 

B 

B 

B 


68 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 

B 101 
B 102 
B 103 
“BIOT 
B 105 
B 106 
B 107 
B 108 
B 109 
B 110 
B 111 
B 112 
B 113 
B 114 
B 115 
B 116- 
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SUBROUTINE EQUAT <DUM1.uUM2,DUM3.5B.EALT,ELAT,ELONG,ERR> c 1 

COMMON FRONT. SHEET 1* SHEET 2 *FSHEET C 2 

COMMON B<2UO),VNl<200>.VN2<200>,VV3<200),ARC(20n>,VNEAR(3).VBQ(3)» C 3 

1VSaVE(3).B0,8NEaR. JUP.MMM C 4 

COMMON VCONJ(3),ALTO C 5 

c c 6 

c SUBROUTINE TRACES field line from a given point to MINIMUM B C 7 

c C 8 

DIMENSION V ( 3 » 3 ) » VN ( 3 ) # VP<3>* R1 ( 3 > > R2<3), R3<3) C 9 

ERR 1 =ERR c 10 

MMM *1 c 11 

jU® = l c 12 

V ( 1 , 2 > = DUMl c 13 

' V<2.2>=DUM2 C 14 

V<3.2>=DUM3 C 15 

ARCIlIsO. C 16 

DCLT = 1 .5708 c 17 

1 AHC ( 2 ) * V ( 1 » 2 ) *SQRT < ERR ) * 0 • 3 C 18 

IF ( V ( 2 » 2 ) - DCLT ) 2.3,3 C_ 19 

2 ' ARC ( 2 1 = * ARC ( 2 ) C 20 

3 CAUL START < R1 . R2 , R3 . B , ARC . ERR • V ) C 21 

IF (JUP.LT.O) GO TO 5 C 22 

DO 4 1=1,3 C 23 

VP ( I ) s V ( I , 2 ) c 24 

4 VN ( I ) = V ( 1 , 3 ) C 25 

CALL LINES <R1.R2.R3,B,ARC,ERR#U.VP.VN) C 26 

IF (J.LT.200) GO TO 7 C 27 

ERR=4 . *ERR c 28 

GO TO 6 C 29 

5 JUP=1 C 30 

gRRiERR*0.1 C 31 

6 CONTINUE ‘ ~C '32 

WRITE (6.8) ERR C 33 

GO TO 1 C 34 

7 eRR*ERR1 C 35 

EB=8NEAR C 36 

EALT=VNEAR(1) C 37 

EIaTs90.-VNEAR<2>*57. 2957795 " ~ — C ' 38 

EL0NG=VNEAR(3)*57. 2957795 C 39 

IF (ELONG.GT.180- > ELONG*ELONG-360 . C 40 

IF (EL0NG.LT. (-180. )) ELONG = ELONG*360 . C 41 

RETURN • C 42 

C_ C 43 

C C ' 44 

C C 45 

8 FORMAT (HX.24H ERROR CHANGED IN EQUAT , E15 . 4 > C 46 

END C 47 
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c 

c. 

c 

c 


1 

2 

3 


4 


5 


C 

C 

6 


SUBROUTINE INSECT (Pl»P2,P3,CALT * CL AT » CL 0 N G # E R R_>_ 

COMMON 8(200 )!vn!C200 ) » VN2<200> #VN3<200)* ARC(200).VNEAR(3)#VBO(3) # 
1VSAVE 13) •BQ*BNEAP*JUP#MMM 
COMMON VC0NJC3) , ALTO 

SUBROUTINE TRACES FIELD LINE D 0 HNW AR D S UNTI L 

ALTO (*1. RE) is REACHED 


DIMENSION V ( 3* 3) * VN<3># VP<3># Rl<3), R2<3), R3(3) 

erri^srr 

J*0 

s 1 T*aBS < SI N [92 )J 

MMM*2 
Vd#2)=Pl 
V ( 2* 2 ) *P2 
V ( 3# 2 ) *P3 
ARC( 1 ) =0 • 

DCLILU5708 _ ... K — 

ARC ( 2 ) S V( i , 2 ) *SqR^ (ERR )*0 *3 
I F ( V ( 2# 2 ) -DCLT ) 3,3,2 
ARC< 2 ) = * ARC ( 2 ) 

CALL START ( Rl, R2 > R3,8, ARCtERR# V ) 

DO 4 1*1,3 

vp(i)*vu,2) _ 

VN(1)«V(I,3) 

call LINES <R 1,R2,R3,B, ARCtERR# J,VP*VN> ___ 

IF (J.LT.200) 60 TO 5 

ERR*4.*ERR 

WRITE (6,6) ERR 

GO TO 1 

“errse'rri 

JUP = J 

CAlT=VC 0NJ<1) 

CLAT=90.-VCONJ(2>*57.2957795 
CLONG*VCONJ< 3) *57. 2957795 

IF (CLONG.GT.10O* >_CLJ^NG*CLOJYG060_._ _ 

IF ( CL0N6 . LTTT*iB OTT ) CLONGsCLONQ^ 360. 

RETURN 


FORMAT (01Xf23HERROR CHANGED IN I NSECT , E16 • 4 ) 
END 


D 1 

D 2 
D 3 
D 4 
D 5 
D 6 
D 7 
D 8 
D 9 
D 10 
D 11 
D 12 
D 13 
D 14 
D 15 
D 16 
D 17_ 
■ Die 

D 19 

D 20 

D 21 

~ D "22 

D 23 

“ D“ 24 
D 25 
D 26“ 
D 27 

D" 20 
D 29 
“ D 30 
D 31 
D 32 
D 33 
D 34 
D 35 
D 36 
D 37 
D 38 
D 39 
D 40 
D 41_ 
D 42 
D 43- 
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(jonon j -4 » ^ » CH M 0 0,0 0 


SUBROUTINE BESECT < RUMl , RUM2 , RUM3 . B ALT , BL AT , 8L0NG » ERR > 

COMMON FRONT, SHEET 1 . SHEET2, FShEET 

COMMON B(200) , VNl < 20 0 > , VN2 ( 20 0 > ,VN3(200) » ARC < 20 0 ) , VNE AR ( 3 > . VBO ( 3 > , 
lVSAVt<3).B0.BNEAR» JUP,MMM 
COMMON VC0NJ(3) , ALTO 

SUBROUTINE traces FIELD line DOWNWARDS FROM A GI_VJ^NJ’0JNT_UNT1L 

A PREFIXED B-VALUE IS REACHED 

DIMENSION V ( 3 » 3 ) » VN(3), VP<3). Rl(3). R2(3>, R3(3) 

ERRlsERR 

MMMS3 

JUP=1 . , 

V < 1 » 2 ) sRUMl 
V(?>2>*RU M2 

V t 3 » 2 ) S RUM3 
A«C(1>*0. 

DCUTsl.5708 

ARC(2)«V<1,2)*SQRT(ERR)*0.3 

IF ( V ( 2 » 2 ) * DCLT ) 3,3,2 
ARC ( 2 > * • ARC ( 2 > 

CALL START ( R1 , R2 . R3 . B , ARC # ERR. V ) 

IF (jUP.LT.O) GO TO 5 _ ' 

DO 4 1*1,3 

VP(I>*V(I,2) _ 

VN< I>*V( 1,3) 

call LINES (R1,R2.R3,B, ARC. ERR# J, VP,VN> 

IF (J.LT.200) GO TO 7 
ERR* 4 • *SRR 
GO TO 6 

JUP*1 _ 

err*err*6. i 

CONTINUE 
write (6,8) ERR 
GO TO 1 

err*erri 

JUPsJ 

BALtrVBOU) 

BLAT* 90. -VBO (2 >*57. 2957795 
BL0NG*VB0(3)*57. 2957795 
IF (BLONG.GT.180. > BLONG*BLONG-360 . 

IF (BLONG.LT. ( -180. > ) BLONG*BLONG*360 . 

RETURN 


FORMAT (#1X,24H ERROR CHANGED IN BESECT, E15. 4) 
END 
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SUBROUTINE INVAR ( DUM1 , DUM2 , DUM3 , ERR , BB, F l ) 


F 

1 

•' • — 

COMMON FRONT, SHEET l . SHEETZ . FSHEET 


F 

2 


COMMON B(200),VNl(200>.VN2<200>» VN3<200>» ARC(200).VNEAR<3) » VB0(3). 

F 

3 


1VSaVE<3)»B0,BNEAR» JUP.MMM 


F 

4 


COMMON VCONJI3), AUTO 


F 

5 

c 



F 

6 

c 

SUBROUTINES INVAR. START, LINES and integ are based on 

MCILWAINS 

F 

7 

c “ 

INVAR coot, "CONVENIENTLY IMPLEMENTED 


F 

8 

c 



F 

9 


DIMENSION V ( 3 , 3 ) , VN(3), VP ( 3 ) » BEG(200 )» BEND<200). 

BLOG ( 200), EC 

F 

10 


10(200), RK3), R2(3). R3(3) 


F 

11 


mmm»o 


F 

12 


JUP*1 . 


F 

13 


V ( 1 , 2 ) =DUM1 


F 

14 


V ( 2, 2 ) S DUM2 


F 

15 


V( 3, 2 > =DUM3 


F 

16 


erri=err 


F 

17 


ARC ( 1 > *0 . 


F 

18 

1 

ARC(2)»V(1,2)*SORT(ERR)*0.3 


F 

19 


DClT *1 . 5708 


r 

20 


IF ( V ( 2 , 2 > • DCLT ) 2,3,3 


F 

21 

2 

ARC ( 2 ) = ■ ARC ( 2 ) 


F 

22 

3 

CALL START (R1,R2,R3,B,ARC« ERR# V > 


F 

23 


IF (JUP.LT.O) GO TO 8 


F 

24 


DO 4 1=1,3 


F 

25 


VP< I ) *V< I , 2) 


F 

26 

4 

VN< I ) S V ( I # 3 ) 


F 

27 


CALL LINES <R1.R2,R3,B,ARC»ERR* J.VP.VN) 


F 

28 


IF (J.LT.200) GO TO 5 


F 

29 


ERR=ERR*4. 


F 

30 


write (6.9) ERR 


F 

31 


GO To 1 


T“ 

32 

5 

err=erri 


F 

33 


JUP = J 


F 

34 


DO 6 Jsl.JUP 


F 

35 


ARC(J)=ABS( ARC(J) ) 


F 

36 

6 

BLOG( J)*ALOG(B< J) ) 


F 

37 


JEP=JUP-1 


F 38 


DO 7 J=2. JEP 


F 

39 


ASUM=ARC( J)*ARC( J*l> 


F 

40 


DX=BLOG(J-l)-BLOG( J) 


F 

41 


DN=ASUM*ARC< J)*aRC( J* l) 


F 

42 


BCO=( (BLOG( J-l)-BUOG( J*l> >*ARC< J)**2-DX*ASUM**2)/DN 


F 

43 


CCO* ( DX* ARC ( J*1 ) * ( BLOG ( J ) 'BLOG < J + l ) ) * ARC f J ) )7 DN 


r 

44 


SA=.75*ARC< J) 


F 

45 


SC=SA* . 25* ASUM 


F 

46 


DCO=BLOG< J-l)-CCO*SA*SC 


F 

47 


ECO( J)=BCO*CCO*<SA*SC) 


F 

48 


BEG( J)=EXP(DCO+ECO( J)*.5*ARC( J) ) 


F 

49 

T” 

"bEnT5('j )=EXP ( DCO*ECO( J)* .5*1 ASUM+ARC < J) i ) 


’ f 5"D 


BEG ( JUP )* BEND ( JEP ) 


F 

51 


BEND( JUP)=8( JUP) 


F 

52 


ECO( JUP)»(2.0/ARC( JUP) )*ALOG<BEND( JUP) /BEG ( JUP) ) 


F 

53 


CALL INTEG (ARC, BEG, BEND, B.JEPiECO, FLINT) 


F 

54 


FlsFLINT 


F 

55 


BB = B ( 2 ) 


F 

56 

8 

CONTINUE 


F 

57 


RETURN 


F 

58 

C 



F 

59 

C 



F 

60 

c 



F 

61 

9 

FORMAT (79X.26H ERROR INCREASED IN INVAR.E15.4) " 


T 

62 


END 


F 

63- 
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SUBROUTINE START (R1#R2»R3«B«ARC,ERR#V) 

COMMON FRONT, SHEET i,SHI=ET2*FSHEET " 

COMMON B(200)#VNl<200),VN2(200>#VN3(200)#ARC(200)#VNEAR(3), VB0(3)# 
1VSAVE(3),0O»BNEAR. JUPiMMM 
COMMON VC0NJ(3),ALT0 
DIMENSION V ( 3 # 3 ) # R1 ( 3 ) # R2<3># R3(3) 

LOOP = l 

SIT=A0S(SIN( V(2,2) ) ) 

IF < V ( 3 # 2 ) ) 2^3,3 
V(3.2>=V(3. 2)^6. 283105307 
GO TO 1 

CALL MODMAG (V<1#2),SIT#V(3, 2)#BR,BT,BP#B(2),V(2,2)) 

R2 ( 1 ) s BR/9 ( 2 ) 
t)N = B < 2 ) *V C l , 2T 
R2 ( 2 ) 5 BT/DN 
R2(3)«BP/(DN*SIT) 


G 1 
' <T “ ~2 


G 

G 

G 

G 

G 

' G‘ 

G 

G 

G 

G 

G 

TT~ 

G 

G 


3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14" 

15 

16 



lS = 0 

G 

17 

4 

1)0 5 1*1.3 

G 

18 

5 

V<1.1)«V(I.2)-ARC(2)*R2(I) 

G 

19 


S I T* aBS < 8 1 N fVT 2 » 1 > ) ) 

■ - ir~Tir 

6 

CALL MODMAG < V . S l T , V < 3, 1 > » BR, BT , BP. B< 1 > , V ( 2. 1 > > 

G 

21 


IF (MMM.LT. 2) GO TO 7 

G 

22 


IF (MMM.eO.2) GO TO 8 

G 

23 


IF (B(l>-8(2>> 10,10.9 

G 

24 

7 

IF (B<l)-B(2>> 9.10.10 

G 

25 

8 

iFTva;i)-V(1.2)) 9,‘971(f ' ■ “ 

G 

”2* 

9 

ARC( 2 > *• ARC ( 2 ) 

G 

27 


IF (LOOP. EG. 2) GO TO 14 

G 

20 


L00P=2 

G 

29 


GO TO 4 

G 

30 

10 

RK1)*BA/B(1> 

G 

31 


ARC( 3 > * ARC ( 2 ) 

G 

”3T 


DNs8 ( 1 ) *V ( 1 . 1 ) 

G 

33 


R1(2)*BT/DN 

G 

34 


R1(3)*BP/(DN*SIT> 

G 

35 


DO U HI, 3 

G 

36 

11 

V( I,l)*V( I,2>-4RC(2)*(R1( I )*R2< 1 >>/2. 

G 

37 


SIT5ABS(SIN(Vt2,IT77 " ~ 

G“ 

“35 


IS*IS*1 

G 

39 


GO TO (6,12), IS 

G 

-40 

12 

DO 13 1*1.3 

G 

41 

13 

V< I.3)*V( 1 , 2 > * ARC t 3 > * ( (1.5)*R2(I )-.5*Rl( I >) 

G 

42 


GO TO 15 

G 

43 

14“ 

JUP*-1 ' — 

G 44 

15 

CONTINUE 

G 

45 


RETURN 

G 

46 


END 

G 

47 


24 



1 

T” 


3 


4 


5 


6 

“7~ 

8 

9 

10 
11 

“12 

13 

14 

15 

T5~~ 
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SUBROUTINE LINES (R1,R2,R3,B# ARC,ERR, j.vp.vn) 

COMMON FRONT, SHEET1,SHEET2#FSHEET 

COMMON B<200 > •VN1<200)*VN2(200> tVN3(200) * ARC(203> » VNEARC3) *VBO(3> # 
1VS AVE < 3 ) .BQ.BNEaR.JUP.MMM 
COMMON vconj<3) , alto 
INTEGER FLAG1.FLAG2 

DIMENSION Rl<3>, R2<3), R3C3), Vn(3>, VP<3). RA(3) 

m*mmm 

flagi=o 

DEL=0.01 
CRE = 0 * 25 

IF CERR-0 .15625) 1.2,2 

CRE=(ERR**0. 333333333) _ __ __ 

A3= ARC ( 3 ) 

A AB- aBS C A3 ) 

SNA=A3/AAB 

Al=ARC(l) 

A2 = ARC ( 2 ) 

A06=A3*A3/6. 0 _ 

VNi(2)=VP(l) 

VN2(2)=VP<2) 

VN3(2)=VP(3) 

J = 3 
I LP»1 

IS = 1 __ 

GO TO 8 
I S = 1 
JSJ + 1 

A06=A3*A3/6. 0 
ARCJ=A1*A2*A3 
AD=( ASUM*A1)/AA 
80= ASUM/BB 
CD=Al/CC 
DO 7 1=1,3 

DD=Rl( I )/A a-R2< I >/BB+R3< I >/CC 
GO TO (5,6), IS 

RT=RK I )•( AD*Rl( I )-BD*R2( I )^CD*R3( I )-DD*ARCJ)*ARCJ 
RAU )-RKD 
Rl( I )=R2( I ) 

R2< I > = R3 ( I > 

R3( I >=RT 
VP( I )=VN< I ) 

RBAR=(R2( I >*R3< I ) >/2.-DD*A06 
VN( I)=VP( I )*A3*RBAR 
IF ( VN ( 2 ) ) 9.10,10 
VN(2)=-VN<2) 

IF (VN(2>*3. 141592653) 12.12#U 
VN(2)=6.2831053O7-VN<2> 

GO TO 10 

IF ( VN < 3 > ) 13,14,14 
VN<3)*VN<3>+6. 283185307 
GO TO 12 

IF (VN(3)-6. 283185307) 16.16, l5 
VN<3 )sVN<3>-6. 283185307 
GO TO 14 

— ‘GO TO "(17,18). IS 
SIT*aBS<SIN< VN(2) ) > 

PRE1=VN<1) 

PRE2=PRE1*VN(2> 
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H 

H 

H 

H 

H 

H 

M 

H 

H 

H 

H 
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H 

H 

H 

H 

W 

H 

H 

H 

H 

H 

' H 
H 
H 
H 
H 
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1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 


TTf 

H 15 


16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 


"TT 32' 
M 33 


34 

35 

36 

37 


H ~ T58 
H 39 


H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 


40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 " 

57 

58 

59 
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18 


19 


20 

21 


22 

23 


PR E3=PRE1*SIT*VN<3) 

call MODMAG (VN,SIT,VN(3),BR,BT,BP#B(J)#VN(2)) 

R3(1)=BR/B( J) 

DN=B( J>*VN(1> 

R3(2)=BT/DM 

R3(3)=BP/<DN*SIT) 

ASUM=A3^A2 

AA = ASUM*A2 _ . 

BB=A3*A2 
CC=ASUM*A3 
IS = 2 
GO TO 4 

SJT=ABS<SIN(VN<2))) 

IF (VM(1) .GT.8. ) GO TO 19 __ 

B ( J) =B < J ) * ( ( PREl/VN ( 1 ) ) **3 ) 

IF (M.EQ.l) GO TO 23 

QRT=.5*ABS<R3(1) >/( .1*ABS<R3<2)*VN<1) ) ) 

XMABS<VN(l)-PRei>^0RT*ABS<VN(iUVN(2)-PRE2)^ABS(VN<l>*SIT*VN(3)-P 
IRE 3) )/( AAB*ERR*SQRTC1.*QRT*QRT> ) 

GO TO _<23#20,23), ILP 

IF (X-3.3) 23 * 21 > 21 

A3=A3*0«2*<8*0*X)/<0.8+X) _ 

J S J-1 

ILP33 

ASUM=A2+A1 

AA = ASUM*A1 _ _ 

B B = A 2 * A 1 

CC=ASUM*A2 

DO 22 1*1*3 
VN ( I ) * VP ( I ) 

R3< I ) «R2< I ) 

R2( I ) *R1 ( I ) _ 

R 1 (I ) s R A (I ) 

GO TO 37 
VNl ( J) *VN( 1) 

VN2( J>=VN(2) 

VN3( J)*VN<3) 

IF (M.E0.2) GO TO 27 

[F^TR . E073T GO TO 29 ~ 

IF <b< J- l) .GT.B< J) ) GO TO 3l 
IF (M.EQ.l) GO TO 25 


H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

“R 

M 

H 

H 

H 

H 

TT 

H 

H 

H 

H 


60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 
_79_ 

80 

81 

82 

83 

84 
05 


86 

87 

08 

89 

90 

91 

~«r 

93 

94 

95 

96 

97 
W 
99 


H 
H 
H 
H 
H 

“FT 
H 

H 100 


IF (FLAGl.EQ.l) GO TO 31 H 101 

FlAGisl H 102 

DO 24 T* 1.3" ” TO 

24 VNEAR< I ) *VN ( 1 ) H 104 

BNEAR«B(J) H 105 

GO TO 31 H 106 

25 BNEARsB(J-l) H 107 

DO 26 1*1,3 _ _ __ 1108 

26 VNEAR ( I ) *i VP Cl) ~ H 109 

GO TO 39 H 110 

27 IF ( VN ( 1 ) . GT • AL TO ) GO TO 31 H 111 

FAC=(ALTO-VP<l) )/(VN<l)-VP(l)> H 112 

DO 28 1*1,3 H 113 

28 VCONj( I )nVP( I )*(VN< I ) - VP ( I >)*FAC H 114 

aRc( j) = arc( j)*faC “ ' "R liB 

VNK J)sVCONJ(l) H 116 

VN2< J)=VC0NJ(2) H 117 

VN3 ( J ) = VCON J ( 3 ) H 118 
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GO TO 39 

29 IP (8(J).L t .B0) GO TO 31 

nr* <bo-b< J-i> )/(B( J>-BU-i)r 

DO 30 1*1.3 

30 VB0( I)*VPC )*(VN< I )-vp(D)*r*c 
aRC(J)=ARC< J)*FaC 

B< J)*BO 

VNl ( J ) *VB0( 1 > 

VN2( j) *VB0i(2) 

VN3 ( J ) * VBO ( 3 ) 

GO TO 39 

31 IP I J.GE.200) GO TO 39 
Al*A2 

JP (M.NE.O) GO TO 32 
1~f~Tb C3 T -SIDTl 27T2 .36 " 

32 ItP*2 
A2 = A3 

IP (M.EQ.l) GO TO 37 
A3=A3*.2*l8.*X)/< .8*X) 
aMs(2.-R3(2)*VN(1) )*VN(1>*CRE 
FTTaBSTaJ’)'* A HI 34,34,33 

33 A3=SNA*AM 

34 IP <SNA*R3<1>*.5> 35.35.37 

35 aM=-.5*SNA*VN(1)/R3<1) 

IF ( ABS ( A3) * AM) 37,37.36 

36 A3=SNA* AM 

T7 ARC < J*lT* A3" ' 

AA0»a8S( A3) 

GO TO 3 

38 CONTINUE 

39 RETURN 

END 


H 119 
H 120 
W 121 
H 122 
H 123 
H 124 
H 125 
H 126 
W 127 “ 
H 128 
H 129 
H 130 
H 131 
H 132 

“H~173 
H 134 
H 135 
H 136 
H 137 
H 138 
H 1J9 
H 140 
H 141 
H 142 
H 143 
H 144 

' TTT25 
H 146 
H 147 
H 148 
H 149 
H ISO- 
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SUBROUTINE 1NTEG CARC.BEG.BEND.B, JEP.ECO.FI ) I 1 

COMMON FRONT. SHEET1.SHEET2.FSHEET 12 

COMMON B( 200 ) . VNl < 200 ) . VN2<200 > , VN3( 200 ) . ARC(200 ) » VNE*R< 3) . VBO< 3) » I 3 

1VS4VE(3).BO,BMEaR. JUP.MMM I 4 

COMMON VCONJ(3) . ALTO I 5 

DIMENSION BEG(200 ). BEND<200>. ECO < 20 0 > l 6 

1 KKsJgP I 7 

IF (KK-4) 3.2.4 I B 

2 • KK*KK-1 I 9 

3 A»3( KK-1 ) /B< 2 ) 1 10 

X2=B<KK)/B<2) i 11 

X3 = B ( KK*1 ) /R ( 2 ) 1 I 2 

aSUM=ARC<KK)*ARC<KK*1) __ I 13 

DN:ARC(KK)*ARC(KK*1)*ASUM I ’ 14 

B8=(-A*ARC(KK*1)*(ARC(KK)*ASUM)+X2*ASUM**2*X3*ARC(KK)**2)/DN I 15 

C=(A*ARC(KK+l)-X2*ASUM*X3*ARC(KK))/DN I 16 

Fl=1.570796326*a.-A*BB*B8/(4.*C))/SORT(ABS(C)) I 17 

RETURN I 18 

4 TsSORT ( l.-BEND(2)/B<2) ) I I 9 

FI = (2.*T-AL0G<<1.*T»/<i.-T)))/ECO(2> ' " T 20 

IF (B(2)-BEND(KK>) 6.6.5 I 21 

K K = KK ♦ 1 I 22 

T=SQRT( ABS(1 . 0-BEG<KK)/B(2> ) ) 1 23 

FlsFl-(2.*T-ALOGUl.*T)/(l.-T) ) )/ECO«KK) I 24 

KKsKK'l I 25 

DO 15 1*3. KK I' 26 

ARG1*1.-BEND(I)/B(2) I 27 

IF (ARG1) 7.7.8 I 28 

7 TE=l.E-5 I 29 

GO TO 9 I 30 

8 TE=SQRT< ARG1) I 31 

9 aRG1=1.»BEG(I)/B<2) - I 32 

IF ( a«G1 ) 11,11,10 ' I 33 

10 T8 = SQfiT ( ARG1 ) I 34 

GO TO 12 I 35 

11 TB*1 . E*5 I 36 

12 IF ( ABSCECOC I ) )-2.E*5) 13.13,14 I 37 

13 ' Fl sFl * { (TE*TB)*(AfiC( I ) ♦ ARC ( 1*1) ) )/4. " " ' ~ 1 38 

GO TO 15 I 39 

14 Fl=FI*(2.*(TE-TB)-AL0G( C 1 . *TE ) * < 1 . -TB)/( <1. -TE)*(1 .*TB) ) ) )/ECO< I > I 40 

15 CONTINUE I 41 

RETURN I 42 

END I 43 
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SUBROUTINE MODMaG (RR.SINTH.PPMI.BR.BTHETa.BPHI.BB.ThET) J 1 

COMMON FRONT « SHEET 1 , SHEET 2* TSHEET J 2 

J 3 

SUBROUTINE assembles MAGNETIC FIELD from tail, magnetopause AND J 4 

INTERNAL DIPOLE - FOR ENQUIRIES WRITE TO GILBERT MEAD. GODDARD J 5 

SPACE FLIGHT CENTER. GREENBELT MARYLAND 20771 J 6 

J 7 

DTmEn51WGS<7,7T ' -J r 

DIMENSION G ( 7 . 7 ) . C0NST<7.7), P(7.7>, DP<7.7), SP<7), CP(7) J 9 

COSTH*COS(THET) J 10 

SiNPHls-SINIPPHI ) J 11 

COSPHU-COSCPPHI ) J 12 

RO=FRONT J 13 

RT* SHEET! - -j - TT - 

R2=SHEET2 J 15 

BCS*FSHEET J 16 

IF (JFIRST-13) 1.5,1 J 17 

JF I RST *13 J 18 

J 19 

SET TJP TNTTTAL CONST*NTS”THE “FIRST TIME” WOUND 3 — 2TT 

J 21 

DO 2 Nsl,7 J 22 

DO 2 M s 1 , 7 J 23 

GG(N,M)«o. " J 24 

NMAX = 7 J 25 

THE FOLLOWING COEFFICIENTS ARE SCHMIDT-NCRMAL IZED J 27 

j 

GG<2,l)*-0.2511lE5 J 29 

GG(3,2)*-0. 12424 E5 J 30 

GG(4,1)»-0.00716E5 J 31 

GGT4,3T«-0. 02333E5 ' J~ 32 

GG(5,2)»-0.02397E5 J 33 

GG(5,4)«-0 . 00163E5 J 34 

GG(6,1)s0.C0569E5 J 35 

GG(6,3)«-0.01078E5 J 36 

GG(6,5)«-0.00103E5 J 37 

GG ( 7 . 2 ) » U . 0 0126E5 J — 3B“ 

GG<7, 4)i-o. 001B7E5 J 39 

GG(7,6)»-0.0004iE5 J 40 

P<1»1>»1.0 J 41 

DP ( 1 , 1 ) » 0 . J 42 

SP<1)*0. J 43 

CP{l)ai. J— 4-4- 

DO 3 N=3.NMAX J 45 

FNsN J 46 

N2 = N-2 J 47 

DO 3 M*1 , N2 J 48 

FM=M J 49 

”3 CONST (M.M) = ( (FN-2. )**2-CFM-l.)4*2r/rC2.*FN v 3.1*T2T4FT<|--r.T) J~ 51T 

DIMENSION SHM 1 DT ( 7 , 7 ) J 51 

SHMlDTd, 1)81.0 J 52 

DO 4 N=2,7 J 53 

FNsN J 54 

SHMIdT(N,1)=SHMIDT(N-1,1)*(FN4.FN-3.0>/(FN-1.0) J 55 

rACT = 2>0 - — — j- 56 

DO 4 Ma2» N J 57 

FM=M J 58 

SHMIDT ( N, M ) =SHM IDT ( N , M-l ) *SQRT ( < FN-FM*1 . 0 ) *F ACT/ ( FN+FM-2 . 0 ) ) J 59 
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11 


12 

13 

14 


15 


FACTal. o 

IF (R0-R00LD) 6,9,6 
R00LD=R0 

DIMENSION F AC ( 7 ) 

FAC<2>*R0**3 
DO 7 Ns 3 , NM AX 
FAC(N>*R0*FAC(N-1> 

DO 8 N*2,NMAX 
DO 8 M«i,n 

G<N, M)=SHMI [)T(N,M)*GG(N,M)/FAC(N> 

CONTINUE 

BEGIN CALCULATION FOR SPECIFIED INPUT 

CTsCOSTH 

STsSINTH 

SP(2)=SINPHI 

CH(2) s C0SPHI 

CALCULATE SIN(M*PHI) AND C0S(M*PHI) 

DO 10 Me3 , NMAX 

SP(M)=SP(2)*CP(M-1 )*CP(2)*SP(M-1) 
CP(M)=CP(2)*CP('M-l)-SP(2)^Sp{M-l) 

R = RR 
A 2 1 . 

RA=R/A 
RO A s 1 . 

BR = 0 . 

BTHETAsQ . 

BPHI = 0 . 

F N s 1 ■ 

CALCULATE SPHERICAL HARMONICS FOR CAVITY FIELD 

DO 16 Na2,NMAX 
SUMRsO . 

SUMT = G . 

SUMP s C . 

FMsfl . 


J 60 
J 61 
* "J 62 
J 63 
J 64 
J 65 
J 66 
J 67 
J 68 
J 69 
J 70 
J 71 
J 72 
J 73 
J 74 
J 75 
J 76 
J 77 
J 78 
J 79 
J 80 
J 81 
J 82 
J 83 
J 84 
J 85 
J 86 
J 87 
J 88 
J 89 
J 90 
J 91 
' J 92 
J 93 
J 94 
J 95 
J 96 
J 97 


J 99 

DEVELOP LEGENDRE FUNCTIONS AND THEIR DERIVATIVES BY RECURSION FORM J 100 


J 101 

DO 15 M s 1 , N J m2 

IF (N-M-l) 13,12^11 J 103 

P<N.M)=CT*P<N-1,M)-C0NST<N/M)*P(N-2>M) "J 104 

DP(N,M) aCT*DP(N-l, M)-ST*P(N~1, M) -CONST < N i M > * DP < N- 2 , M ) J 105 

GO TO 14 J 106 

P(N,M)rCT*P(N-l,M) j 107 

DP(N,M)*CT*DP(N-l,M)-ST*P(N-l,M) J108 

GO TO 14 j 109 

P(NiN>=ST*P<N-l,N-l> ~ J HQ 

DP(N,N)«ST*DP<N-1.N-1)^CT*P(N-1.N-1) j hi 

CONTINUE J U2 

TS = G(N,M)*CP(M) j h3 

SUMRsSUMR*P(N,M)*TS J 114 

SUMT=SUMT+DPCN,M)*TS j n* 

SUMP=SUMP^FM*P(N.M)*G(N,M)*SP(H) ' - j 116 

FMzFM^l. j 117 

CONTINUE J US 

BR=9R-ROA*FN*SUMR j 119 
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J 120 
J 121 
J 122 
J 123 
J 124 
J 125 
J 126 


CALCULATE tail FIELD J 127 

~ J128 

RCT*R*CT J 129 

RCT2=RCT**2 J 130 

RSC*R*ST#CP(2) J 131 

T0P*R2*RSC J 132 

B0T=R1*RSC J 133 

Hr^lCS*“(ATANTTT)P/ITCnT-ATAN^0T/RCn>/3.14159265 J134 

0 PH I =BPHI ♦BX*SP( 2 ) . J 135 

BRH0='BX*CP(2> J 136 

BY=BCS*ALOG< <RCT2*T0P**2>/<RCT2*B0T**2> 1/6.28318531 J 137 

BR=BR*BRHO*ST-BY*CT J 138 

BTHET4sBTHETA*8RH0*CT*BY*ST J 139 

" ' J I4T 

ADD DIPOLE FIELD TO CAVITY FIELD j 141 


J 142 
J 143 
J 144 
J 145 
J 146' 
J 147 
J 148 
J 149 
J 150 
J 151- 


R3sR**3 

BR=BR-62000.*COSTH/R3 

BTHETA*BTHETA-31000.*SINTH/R3 

B'R=BR*TVE"* g"5 — 

BTHETAsBTHETA*l.E-05 
BPH I =BPHI *1 • E-05 

BBsSORT<BR*BR*BTHETA*BTHETA*BPHI*BPHI ) 

RETURN 

END 


BTHETA=BTHETA-ROA*SUMT 
BPHI=BPHI*ROA*SUMP 
UDAsROA*R A 

fn=fn*i. 

CONTINUE 
BPH I =8PH I /ST 


